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1.  Introduction 


Sonic  fatigue  has  been  considered  as  one  of  the  major  design  considerations  for  the  Joint 
Strike  Fighter  (JSF).  In  addition,  the  surface  panels  of  many  high-speed  flight  vehicles 
(e.g.,  the  X-33,  RLV,  X-38,  and  X-43,  etc.)  presently  under  development  will  be  exposed 
to  high  levels  of  acoustic  pressure  and  elevated  temperatures.  This  brings  an  urgent  need 
for  the  sonic  fatigue  analysis  and  design  methods  for  aircraft  and  spacecraft  structural 
panels. 

For  safety  and  reliability,  the  design  of  modem  structures  such  as  skyscraper  buildings, 
constructions  housing  nuclear  reactors,  and  naval  and  aerospace  structures  must  take  into 
consideration  various  intense  random  excitations.  These  excitations  include  seismic 
motions  of  earthquakes,  pressure  waves  of  explosions  or  blasts,  jet  noise,  and 
atmospheric  turbulence.  The  three  important  aerospace  problems  of  random  vibration  [1] 
are  the  following:  buffeting  of  aircraft  by  atmospheric  turbulence,  sonic  fatigue  of  aircraft 
and  spacecraft  panels  due  to  jet  noise  impingement  or  boundary  layer  pressure 
fluctuations,  and  the  reliability  of  payloads  in  rocket-propelled  vehicles. 

1.1  Nonlinear  Random  Vibration  Analysis  Techniques 

Stochastically  excited  linear  systems  have  been  studied  in  great  detail  and  numerous 
analytical  techniques  exist  for  both  stationary  and  nonstationary  problems. 
Unfortunately,  the  majority  of  structural  responses  are  nonlinear  and  not  many  techniques 
exist  for  the  analysis.  Crandall  and  Zhu  [1],  To  [2],  Roberts  [3],  and  Spanos  and  Lutes 
[4]  have  presented  excellent  and  comprehensive  reviews  on  techniques  for  nonlinear 
random  vibrations.  The  various  approaches  are  given  briefly  in  the  following  paragraphs. 

1.1.1  Fokker-Planck-Kolmogorov  (FPK)  Equations  Approaches 

The  FPK  equations  approaches  give  an  exact  solution  for  a  restricted  class  of  simple 
problems.  The  most  general  extension  of  FPK  equations  approaches  to  nonlinear  second 
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order  equations  was  due  to  Caughey  [5].  Exact  steady-state  solutions  of  rather  wide  class 
of  Multi-Degrees-of-Freedom  (MDOF)  nonlinear  systems  to  white  noise  are  available 
[6,7].  In  general,  the  transitional  Probability  Density  Function  (PDF)  cannot  be  found 
with  the  FPK  equation  approach.  Without  this  transitional  probability,  it  is  generally 
impossible  to  obtain  the  correlation  function  and  Power  Spectral  Density  (PSD)  of  the 
response.  The  difficulty  in  dealing  exactly  with  solutions  of  stochastically  excited 
nonlinear  systems  has  led  to  an  intensified  effort  to  develop  approximate  methods,  to 
tackle  a  broader  class  of  problems  than  presently  possible  with  the  exact  analysis. 

1.1.2  Perturbation  Approaches 

In  this  approximate  method,  the  stochastically  exited  nonlinear  system  is  treated  in  the 
same  way  as  a  deterministically  excited  system.  The  solution  is  represented  as  an 
expansion  of  the  powers  of  a  small  parameter  which  specifies  the  size  of  the  nonlinearity. 
The  perturbation  approach  was  applied  to  a  continuous  nonlinear  system  by  Lyon  [8]  and 
to  discrete  nonlinear  systems  by  Crandall  [9].  The  perturbation  approximation;  however, 
will  not  give  accurate  result  for  systems  of  large  nonlinearity  [10]  as  shown  in  Figure  1. 

1.1.3  Equivalent  Linearization  (EL)  Approaches 

The  EL  approaches  technique  is  based  on  the  concept  of  replacing  the  nonlinear 
system  by  an  equivalent  linear  system  such  that  the  difference  between  the  two  systems  is 
minimized.  Basically,  the  method  is  the  statistical  extension  of  well-known  Krylov- 
Bogoliubov  equivalent  linearization  method  for  deterministic  vibration  problems.  The 
extension  of  this  approximate  method  to  problems  of  random  excitations  was  made 
independently  by  Booton  [1 1]  and  Caughey  [12].  Atalik  and  Utku  [13]  have  presented  a 
direct  and  generalized  procedure  of  the  equivalent  linearization  approach  for  the  MDOF 
nonlinear  systems  that  may  be  nonlinear  in  inertial,  velocity,  and  restoring  forces.  The 
coefficients  of  the  equivalent  linear  system  can  be  obtained  by  direct  application  of 
partial  differentiation  and  expectation  operators  to  the  functionals  involving  nonlinear 
terms.  For  mathematical  derivations  of  the  equivalent  linearization  technique  and  its 
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applications  to  the  variety  of  nonlinear  dynamic  systems,  the  readers  are  referred  to  the 
published  book  by  Roberts  and  Spanos  [14]. 

1.1.4  Numerical  Simulation  Approaches 

Simulation  or  the  Monte  Carlo  method  estimates  the  response  statistics  of  randomly 
excited  nonlinear  structural  systems  [15-17].  In  the  past,  both  analog  and  digital 
computational  systems  have  been  used  for  Monte  Carlo  simulations.  Only  digital 
systems  are  used  presently.  The  approach  mainly  consists  of  generating  a  large  number 
of  sample  excitations,  calculating  the  corresponding  response  samples,  and  processing  the 
desired  response  statistics.  Obviously,  this  approach  can  be  used  for  estimating  the 
response  statistics  of  both  stationary  and  nonstationary  excitations.  The  major  drawback 
of  this  approach  is  the  computation  time  and  cost. 

The  various  analysis  techniques  discussed  for  nonlinear  random  vibration  systems  in 
Section  1 . 1  were  not  dealt  with  thermal  environment.  A  brief  review  of  sonic  fatigue 
analysis  and  design  methods  for  aircraft  and  spacecraft  structural  panels  in  a  combined 
thermal  acoustic  environment  is  presented. 

1.2  Nonlinear  Random  Response  of  Panels  in  an  Elevated  Thermal 
Environment 

Sonic  fatigue  design  guides  have  been  developed  for  metallic  structures  by  Rudder  and 
Plumblee  [18]  and  for  graphite-epoxy  composite  structures  by  Holehouse  [19].  The 
design  guides  were  based  on  the  semi-empirical  test  data  or  the  simplified  single-mode 
Miles’  approach. 

Vaicaitis  and  his  co workers  have  used  the  Galerkin’s  method  (to  Partial  Differential 
Equations  (PDE)  and  modal  approach)  in  conjunction  with  the  time  domain  Monte  Carlo 
numerical  simulation  [15-17]  for  the  prediction  of  nonlinear  response  of  isotropic  [20,21] 
and  composite  [22,23]  panels  subjected  to  acoustic  and  thermal  loads.  Lee  [24-26]  has 
used  the  PDE/Galerkin  method  in  conjunction  with  the  equivalent  linearization  [14] 
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technique  and  investigated  the  thermal  effects  on  the  dynamics  of  vibrating  isotropic 
plates  and  the  improvement  of  variance  and  cumulants  using  the  abridged  Edgeworth 
series  [27].  The  use  of  the  PDE/Galerkin  method,  however,  limits  its  applicability  to 
simple  panel  planform  of  rectangular  shape  and  simple  boundary  conditions. 

Extension  of  the  Finite  Element  (FE)  method  to  nonlinear  response  of  isotropic  beam 
[28]  and  plate  [29]  structures  under  combined  acoustic  and  thermal  loads  was  first 
reported  by  Locke  and  Mei  using  the  EL  technique  with  an  iterative  scheme.  The 
application  of  the  FE/EL  procedure  was  further  extended  to  composite  panels  by  Mei  and 
Chen  [30].  In  the  FE/EL  solution  procedure,  the  thermal  postbuckling  or  thermal  finite 
deflection  structural  problem  is  solved  first.  The  thermal  deflection  and  thermal  stresses 
are  treated  as  known  preconditions  for  the  subsequent  random  response  analysis.  The 
random  response  thus  considers  only  one  of  the  two  coexisting  thermal  postbuckling 
positions  [31].  The  FE/EL  method,  therefore,  does  not  give  accurate  predictions  for 
snap-through  (or  oil-canning)  and  large-amplitude  nonlinear  random  motions. 
Experiments  by  Ng  and  Clevenson  [32],  Istenes  et  al.  [33],  and  Murphy  et  al.  [34,  35] 
have  shown  that  the  dynamic  response  of  acoustic  excited  thermally  buckled  plates  may 
exhibit  the  following  two  types  of  motion:  (i)  small  amplitude  vibrations  about  one  of  the 
coexisting  static  equilibriums,  and  (ii)  large  amplitude  nonlinear  snap-through 
oscillations  between  and  over  the  two  postbuckling  positions.  Reviews  of  large 
deflection  analysis  in  sonic  fatigue  design  were  given  by  Mei  and  Wolfe  [36],  Benaroya 
and  Rebak  [37],  Vaicaitis  [20,  21],  Clarkson  [38],  and  Wolfe  et  al.  [39], 

This  report  presents  an  efficient  finite  element  method  for  the  prediction  of  nonlinear 
random  response  of  composite  panels  at  elevated  temperatures.  The  system  equations  of 
motion  are  first  derived  in  the  structural  node  degree  of  freedom  (DOF),  then  they  are 
reduced  to  a  set  of  coupled  nonlinear  modal  equations.  Numerical  integration  is  used  to 
obtain  the  panel  response.  The  following  three  types  of  motions  can  be  predicted:  (i) 
linear  random  vibration  about  one  of  the  thermally  buckled  position,  (ii)  snap-through 
between  the  two  buckled  positions,  and  (iii)  nonlinear  random  vibration  over  the  two 
thermally  buckled  positions.  Examples  are  given  for  an  isotropic  plate  and  a  composite 
plate. 
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Figure  1.  RMS  Responses  of  a  Hardening  System  by 
Perturbation,  EL  and  FPK  Approaches  (Iwan  &  Yang,  1972) 


2.  Finite  Element  Formulation 


The  governing  nonlinear  equations  of  motion  are  derived  for  an  arbitrarily  laminated 
composite  plate  subjected  to  a  set  of  simultaneously  applied  thermal  and  acoustic  loads. 
The  thermal  load  is  taken  to  be  an  arbitrary  distribution  and  steady-state,  i.e., 
AT=AT(x,y,z).  The  acoustic  excitation  is  assumed  to  be  a  band-limited  Gaussian  random 
noise  and  uniformly  distributed  over  the  structural  surface. 

2.1  Equations  of  Motion  in  Structural  Node  DOF 

The  element  displacements  are  expressed  in  terms  of  the  node  DOF  as 
u(x,y,t)  =  [_Hu\{wm} 

v(x,7,0  =  L^vJK,}  (1) 

w(x,y,t)  =  [_Hw\{wb} 

where  u,  v,  and  w  are  the  in-plane  and  transverse  displacements  of  the  middle  surface;  the 
vectors  {wm}  and  {w*}  denote  the  in-plane  and  bending  node  DOF;  and  \_HU  J ,  [_//,,  J,  and 
\_HW  J  denote  in-plane  and  transverse  displacement  functions,  respectively. 

The  finite  element  formulation  is  based  on  the  von  Karman  large  deflection  and  the 
laminated  plate  theories.  The  strain-displacement  relations  are  given 
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where  {f0}  and  {k}  denote  the  in-plane  strain  and  curvature  strain  vectors;  and  [9]  is  the 
slope  matrix,  respectively.  The  matrices  [Bm],  [Be],  and  [Bb\  are  the  strain  interpolation 
matrices  corresponding  to  in-plane,  large  deflection,  and  bending  strain  components, 
respectively.  The  subscripts  m,  9,  and  b  denote  that  the  strain  components  are  due  to 
membrane,  large  deflection,  and  bending,  respectively;  and  the  comma  denotes  the 
derivative. 

The  linear  constitutive  relations  for  the  kth  orthotropic  lamina  in  the  principal  material 
coordinates  are 
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where  [ Q ]  is  the  reduced  stiffness  matrix  of  the  composite  lamina,  and  {a}k  is  the 
coefficient  of  thermal  expansion.  The  terms  in  [0]  can  be  evaluated  as  follows: 
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Consider  the  composite  lamina  having  an  arbitrary  orientation  angle  <f>,  the  stress  and 
strain  transformation  relations  from  the  principal  directions  1,  2  to  x,  y  body  directions 
are: 
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where 

[TM)]=  s2  c2  -2sc  I,  fcWH  s’  C1  -sc  I  (9) 

with  c=cos((j)),  .y=sin(<f>). 

Thus,  the  stress-strain  relations  for  a  general  kth  lamina  with  an  orientation  angle  <|>  and  a 
temperature  change  become: 
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or 


w* =bi 

where  [q  ^  the  transformed  reduced  stiffness  matrix  is  given  by, 

[el =[r.(«r'iein(«] 

The  resultant  forces  and  moments  per  unit  length  are 


(11) 


(12) 
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»  [M})=  JW*ou)flk 


and  the  constitutive  equations  for  a  laminate  can  be  written  as 


N)  \  A  B\\s°  \N 


M  \B  D 


where  [A],  [5],  and  [D]  are  the  laminate  extensional,  extension-bending,  and  bending 
stiffness  matrices,  respectively.  The  vectors  {Nat}  and  {Mat}  are  the  in-plane  force  and 
moment  due  to  thermal  expansion 


({iV„}  ,  ‘jtc],ArH(l,z>fe  (15) 

-hi  2 

Using  the  principle  of  virtual  work,  the  element  nonlinear  equations  of  motion  are 
derived  with  the  internal  and  external  virtual  work  as 

A  (16) 
5Wexl  =  J  [5w[p{x,  y,  t )  -  phw„ )  -  8u(phu  „  )  -  Sv(phv „  )}<£4 

A 

and  the  element  equations  of  motion  can  be  expressed  as 
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The  element  matrices  and  load  vectors  are  listed  in  Appendix  A.  The  subscripts  B,  NAT, 
Nm  and  NB  denote  that  the  corresponding  stiffness  matrix  is  due  to  the  laminate 

extension-bending  [5],  in-plane  force  components  {NAT},  {Nin }  =  [A]{zl } ,  and 

{Nb }  =  [B]{k} ,  respectively. 

Assembling  all  the  elements  and  taking  into  account  the  kinematic  boundary  conditions, 
the  system  equations  of  motion  in  structural  node  DOF  can  be  expressed  as 


0  M 
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or 


(18a) 


[M]lfr}+  ([*]  -  ]  +  [Jf,  ]  +  [K2  ]){W }  =  {/>„  }+{;>(/)}  (18b) 

where  [M\,  [X],  and  {P}  denote  the  system  mass,  linear  stiffness  matrices,  and  load 
vector,  respectively;  and  [. K{\  and  denote  the  first-order  and  second-order  nonlinear 
stiffness  matrices  which  depend  linearly  and  quadratically  on  displacement  {W}. 

For  a  given  temperature  rise  AT,  Equation  (18)  can  be  solved  by  numerical  integration  in 
the  structural  node  DOF  with  simulated  random  loads.  This  approach  has  been  carried 
out  for  random  response  analysis  with  simulated  random  loads  by  Green  and  Killey  [40] 
and  Robinson  [41].  It  turned  out  to  be  computationally  costly  because  of  the  following: 
(i)  the  large  number  of  DOF  of  the  system,  (ii)  the  nonlinear  stiffness  matrices  [Ki]  and 
\Ko\  have  to  be  assembled  and  updated  from  the  element  nonlinear  stiffness  matrices  at 
each  time  step,  and  (iii)  the  time  step  of  integration  has  to  be  extremely  small. 
Consequently,  Equation  (18)  is  transformed  into  a  set  of  truncated  modal  coordinates 
with  rather  small  DOF. 
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2.2  Equations  of  Motion  in  Modal  Coordinates 


2.2.1  Symmetric  Laminates 

For  symmetrically  laminated  composite  and  isotropic  panels,  the  laminate  coupling 
stiffness  [5]  is  null  and  the  two  submatrices  in  Equation  (18)  are 

[*d =[*,*.]=«  (19> 

By  neglecting  the  membrane  inertia  term,  the  membrane  displacement  vector  can  be 
expressed  in  terms  of  the  bending  displacement  vector  as 


{*'„  }  =  [K.  ]■'  ({/Ur  }-  K* M  !)  (20) 

Then  Equation  (18)  can  be  written  in  terms  of  the  bending  displacement  as 

[M„  P„ }+  ([JC,  ]  -  [KmT  DR } + [*,„  R  r  Rr } 

(21) 

+ Rrl+  {n(0> 

In  the  above  system,  the  nonlinear  stiffness  matrices  can  be  expressed  in  terms  of  modal 
coordinates.  This  is  accomplished  by  expressing  the  panel  response  as  a  linear 
combination  of  some  base  functions  (modal  transformation)  as 


W=l>,(0fcrMA} 

r-\ 


(22) 


where  corresponds  to  the  normal  modes  of  the  linear  vibration  problem 


0>2r[MMr)=[Kb]ti>br 


(23) 
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The  nonlinear  stiffness  matrices  [ K/bm ]  and  [AT?*]  are  both  in  function  of  {IfV}.  They  can 
be  expressed  as  the  sum  of  products  of  modal  coordinates  and  nonlinear  modal  stiffnesses 
matrices  as 

r=l 

<24) 

r~\  s= 1 

where  the  super  indexes  of  those  non-linear  modal  stiffness  matrices  denote  that  they  are 
assembled  from  the  corresponding  element  non-linear  stiffness  matrices  (see  Appendix 
A).  Those  non-linear  stiffness  matrices  are  evaluated  with  the  corresponding  element 
components  {wb}(r)  obtained  from  the  known  system  mode  {^}(r). 

The  nonlinear  stiffness  matrix  [/f/vm]  is  linearly  dependent  on  the  displacement 
{Wm}.  Recalling  the  membrane  displacement  vector  of  Equation  (20) 


\  r= 1  5=1 

=  t)qs(oy>j 


(S) 


\(rs) 


r=l  5=1 


(25) 


It  is  observed  that  [Ki^m]  is  the  sum  of  two  matrices,  the  first  [A^vmlar is  evaluated  with 
{Wm}AT  (=[Am]-'  {PmAT})  and  the  second  [K2Nm\  is  evaluated  with  {</>J(rs) HKm]'] 
l Kimb](r)m(s) )  as 


[KlNm] = [knX  Qs(0[K2NMr] 

/--•l  5=1 


(26) 


12 


Introducing  a  structural  modal  damping  2 i;rcorM r[I],  where  the  modal  damping,  £r,  can 
be  determined  experimentally  or  pre-selected  from  a  similar  structure.  The  equations  of 
motion  Equation  (21)  are  reduced  to  a  set  of  coupled  modal  equations  as 


[m]  {q}+2Zra,rMrin  fei+fcl+M  M={p} 


(27) 


where  the  diagonal  modal  mass  is 


[m]=W[m(]  w =M,m 


(28) 


the  linear  and  cubic  terms  are 


MM- M  M+OTIaK.  (29) 

/■=! 


[*»- wei  M.k.r  w  {?}  o»> 

r=1  s=l 


and  the  modal  thermal  and  random  load  vector  is 

{?}=Mr  ({/>„} +{/>»(/)})  (3i) 

The  nonlinear  random  response  for  a  given  symmetric  composite  or  isotropic  panel  at 
certain  temperature  can  thus  be  determined  from  Equation  (27)  by  numerical  integration 
scheme.  The  advantages  of  using  the  modal  equations  are  the  following:  (i)  the  number 
of  modal  equations  is  small  (DOF  of  {#}«DOF  of  {Wb}),  (ii)  there  is  no  need  to 
assemble  and  update  the  nonlinear  cubic  term  at  each  time  step  since  all  the  nonlinear 
modal  stiffness  matrices  are  constant  matrices,  and  (iii)  the  time  step  of  integration  could 
be  larger. 
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2.2.2  Unsymmetric  Laminates 


For  unsymmetrically  laminated  composite  panels,  the  laminate  coupling  stiffness  [B]* 0 
which  leads  to  the  two  submatrices  [Kb\  and  [Kikb]  both  are  not  zero.  The  linear 
vibration  from  Equation  (18)  becomes 


t t) 


0 


~k„  ^i| a  r 
a  ki 


(32) 


The  bending  {^}<r)  and  in-plane  modes  are  thus  coupled. 

Follow  the  procedures  for  the  symmetric  laminates,  the  panel  response  is  expressed  as 

The  nonlinear  stiffness  matrices  [ATy]  and  [Ki]  can  be  expressed  as  the  sum  of  products  of 
modal  coordinates  and  nonlinear  modal  stiffness  matrices  as 


[*,]-£»,(«> 


r=l 


\K\Nm 


<*.r 

[«,- 


[*, 


1  bm 


(A>r 

0 


and 


r= 1  5=1 


(34) 


(35) 


where  the  super  indexes  of  those  nonlinear  modal  stiffness  matrices  denote  that  they  are 
assembled  from  the  corresponding  element  nonlinear  stiffness  matrices.  The  element 
nonlinear  stiffness  matrices  are  evaluated  with  the  corresponding  element  components 
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{w6}(r)  and  {Hv}(r)  obtained  from  the  known  system  modes  {^}(r)  and 
respectively. 

With  the  introducing  a  structural  modal  damping  2 <%ro)rMr[I],  equations  of  motion 
Equation  (18)  reduce  to  a  set  of  coupled  modal  equations  as 


M  {ij}+24rarMrii)  fe}+(fc]+k]+M  w=m 


(36) 


where  the  diagonal  modal  mass  and  linear  stiffness  matrices  are 


(M  fc  ]) =  Wr  (M  [x  -  xmT  Ms*] 


(37) 


the  quadratic  and  cubic  terms  are 


r=l 

(38) 

M?}= w±±™\K»rm 

(39) 

r= 1  5=1 


and  the  modal  thermal  and  random  load  vector  is 

i?!=Mr  (ta-}+ {p(o9  (4°) 

The  nonlinear  random  response  for  a  given  unsymmetric  composite  panel  at  certain 
temperature  can  then  be  determined  from  Equation  (36),  using  numerical  integration. 
The  advantages  are  as  follows:  (i)  DOF  of  {q}  is  small,  (ii)  no  need  to  assemble  and 
update  the  nonlinear  quadratic  and  cubic  terms,  and  (iii)  the  time  step  of  integration  could 
be  larger. 
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2.3  Stress  and  Strain  Calculations 


After  the  modal  displacement  {q}  for  a  given  combination  of  acoustic  load  and  elevated 
temperature  case  is  determined,  {Wb}  and  {Wm}  can  be  evaluated  with  Equation  (22)  and 
Equation  (25)  for  the  symmetric  and  Equation  (33)  for  the  unsymmetric  panels, 
respectively.  The  element  in-plane  strain  {s0}  and  curvature  {k}  can  be  calculated  using 
Equation  (3)  to  Equation  (5),  respectively.  The  element  strains  are  then  obtained  from 
Equation  (2),  and  stresses  for  the  kth  layer  are  obtained  using  Equation  (10).  Stress  and 
strain  in  the  material  principal  directions  are  then  obtained  using  Equation  (8).  For  the 
displacement  based  finite  element  method,  the  stress/strain  calculation  is  not  as  accurate 
as  displacement  calculation.  According  to  Barlow  [42],  and  Cook  et  al.  [43],  the  stresses 
at  Barlow  points  are  calculated  and  the  result  is  extrapolated  to  the  nodal  points  or  other 
desired  points.  The  global  stresses  and  strains  are  averaged  from  different  local  nodal 
values,  which  share  the  same  global  node  number. 
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3.  Results  and  Discussions 


3.1  Uniform  Distribution  Random  Pressure 

Consider  a  random  pressure  p(x,y,t)  acting  on  the  surface  of  a  high-speed  flight  vehicle. 
The  pressure  acting  normal  to  the  panel  surface  varies  randomly  in  time  and  space  along 
the  surface  coordinates  x  and  y.  The  pressure  p(x,y,t)  is  characterized  by  a  cross-spectral 
density  function  Sp( p,  co),  where  xi  -  x2  and  p=yi  -  y2  are  the  spatial  separations 
and  to  is  the  frequency  in  rad/sec.  The  simplest  form  of  the  cross-spectral  density  is  the 
truncated  Gaussian  white  noise  pressure  uniformly  distributed  with  spatial  coordinates  x 
and  y 


s,(S,n,f)= 


S,  if  0  <f<f„ 

0  if  f  <  0  or  f>fu 


(41) 


where  So  is  constant  and  fu  is  the  upper  cut-off  frequency  in  Hz.  The  expression  for  So 
can  be  written  as  [18,22] 


S0  =  p20  \OSPLno  (42) 

where  po  is  the  reference  pressure,  po  =  2.90075  10’9  psi  (20pPa),  and  sound  pressure 
level  (SPL)  is  expressed  in  decibels  (dB).  A  typical  simulated  random  load  at  90  dB  SPL 
is  shown  in  Figure  2.  The  band-limited  white  noise  is  generated  by  a  Fortran  code  given 
in  Appendix  B  that  simulates  a  random  pressure  using  complex  numbers  with 
independent  random  phase  angles  uniformly  distributed  between  0  and  2n.  The  PSD 
value  of  the  random  process  is  obtained  by  taking  the  ensemble  average  of  the  Fourier 
transform  of  the  random  load.  The  PSD  value  is  then  compared  to  the  exact  one  given  by 
Equation  (42).  The  analyses  presented  are  obtained  for  a  cut-off  frequency  of  1024  Hz 
for  the  isotropic  and  the  composite  plates.  The  selected  frequency  bandwidth  is  Ao>=0 
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rad/sec  (1  cycle/sec)  with  the  random  load  prescribed  in  decibels.  For  instance,  for  a 
uniformly  Gaussian  random  load  of  100  dB  over  a  frequency  range  of  0-1024  Hz 
corresponds  to  an  overall  sound  pressure  level  (OASPL)  of  130  dB. 

3.2  Finite  Element  and  Validation 

The  nonlinear  element  equations  developed  in  Equation  (17)  are  general  in  the  sense  that 
they  are  applicable  for  beam  [28],  rectangular  [29,  31,  and  41],  and  triangular  [30,  44, 
and  45]  plate  finite  elements.  The  finite  element  employed  in  the  present  study  is  the  C1- 
conforming  rectangular  plate  element.  The  linear  stiffness  and  mass  matrices  are 
developed  by  Bogner-Fox-Schmit  (BFS)  [46].  The  thermal  geometrical  matrix,  thermal 
load  vectors,  and  nonlinear  stiffness  matrices  are  generated  using  the  expressions  given  in 
Appendix  A.  The  BFS  element  has  a  total  of  24  DOF,  16  bending  DOF  { vty } .  and  8  in¬ 
plane  DOF  {wm}  as  shown  in  Figure  3. 

Accurate  nonlinear  analytical  multimode  results  and  test  data  for  panels  under 
acoustic  and  thermal  loads  are  not  available  in  the  literature.  Validation  of  the  present 
nonlinear  modal  formulation  will  thus  consists  of  the  following  two  parts:  (i)  nonlinear 
free  vibrations  to  assess  the  accuracy  of  the  left  side  of  Equation  (27)  with  zero  damping, 

and  (ii)  linear  random  vibrations  to  validate  the  simulated  random  load  {/’  j  at  the  right 
side  of  Equation  (27)  with  AT=0.  The  accuracy  of  the  nonlinear  stiffness  matrices  in 
modal  coordinates  has  been  verified  by  Shi  et  al.  [47]  for  nonlinear  free  vibration  of 
fundamental  and  higher  modes  of  plates  and  beams.  The  validation  of  simulated  random 
loads  is  by  comparison  of  the  linear  displacements  with  linear  analytical  results  shown  in 
Table  1.  Linear  analytical  random  response  is  given  in  Appendix  C.  The  FPK  method 
[48]  is  an  exact  solution  [49]  to  the  single  DOF  forced  Duffing  equation.  The  present 
time  domain  numerical  simulation  results  are  also  shown  in  Table  1.  The  natural 
frequencies  of  the  lowest  seven  modes  in  increasing  order  are  given  in  Table  2. 
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Table  1.  Comparison  of  RMS  Wmax/h  for  a  Simply-Supported 
15x12x0.040  in.  Isotropic  Plate 


SPL 

(dB) 

Linear  Analytical 

4  modes  7  modes 

FE/L/NS 

4  modes  Err.% 

FPK  [48, 49] 

1  mode 

FE/NL/NS 

4  modes 

90 

0.2759 

0.2759 

0.2760 

0.0362 

0.249 

0.266 

100 

0.8725 

0.8725 

0.8728 

0.0362 

0.592 

0.578 

110 

2.7590 

2.7590 

2.7600 

0.0362 

1.187 

1.432 

120 

8.7248 

8.7250 

8.7281 

0.0362 

2.200 

2.572 

FE:  Finite  Element;  L:  Linear;  NL:  Non-Linear;  NS:  Numerical  Simulation 


Table  2.  Frequencies  (Hz)  of  a  Simply-Supported  15x12x0.040  in.  Isotropic  Plate 


Mode 

(14) 

(34) 

(1,3) 

(34) 

(54) 

(5,3) 

(1,5) 

Exact 

44.078 

181.68 

259.09 

396.70 

456.90 

671.91 

689.12 

FE 

44.082 

181.70 

259.10 

396.75 

457.02 

672.06 

689.29 

Mesh  size  is  10x10  in  a  quarter  plate. 


3.3  Simply  Supported  Isotropic  Plate 

A  simply  supported  isotropic  plate  with  immovable  in-plane  conditions 
u(0,y)=u(a,y)=v(x,0)=v(x,b)=0  is  studied  in  detail.  The  plate  is  of  14  by  10  by  0.04  in. 
(35.6  by  25.4  by  0.1  cm)  and  is  modeled  with  a  14  by  10  mesh  (140  BFS  elements)  in  a 
quarter  plate.  The  number  of  structural  node  DOF  {Wb}  is  560  for  the  system  equations 
given  in  Equation  (21).  The  material  properties  are  E=10.5E7  psi(73  GPa),  v=0.3,  and 
p=2.588xl0'4  lbf-sec2/in.4  (2763  kg/m3).  A  proportional  damping  ratio  of  £ra>r=^s©s  with 
£,1=0.02  is  used.  The  lowest  seven  natural  frequencies  are  given  in  Table  3. 
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Table  3.  Frequencies  (Hz)  of  a  Simply-Supported  14x10x0.040  in.  Isotropic  Plate 


Mode 

(1,1) 

(3,1) 

(1,3) 

(3,3) 

(5,1) 

(5,3) 

(1,5) 

Exact 

58.116 

215.19 

365.98 

523.05 

529.33 

837.19 

981.70 

FE 

58.116 

215.19 

365.98 

523.06 

529.36 

837.22 

981.94 

The  number  of  modal  coordinates  to  be  included  in  the  analyses  for  converged  deflection 
and  strain  is  studied  first.  The  Root  Mean  Square  (RMS)  maximum  non-dimensional 
deflection  and  the  RMS  maximum  strain  versus  number  of  modes  at  SPL  of  120  dB  using 
1,  2,  4,  6,  and  7  modes  are  shown  in  Figure  4.  The  maximum  strain  is  ey  at  the  plate 
center.  The  results  show  that  four  modes  will  give  converged  deflection  and  strain 
responses. 

Two  other  studies  for  accurate  and  converged  response  predictions  were  also 
performed.  They  are  the  finite  element  mesh  sizes  and  the  integration  time  steps.  For  a 
four-mode  solution,  it  was  found  that  a  quarter  plate  model  of  14  by  10  mesh  size  is  more 
than  adequate.  The  time  step  of  integration  At=l/81 92=1. 2207x1  O'4  sec  was  first 
selected,  then  the  time  step  was  cut  into  one-half  with  At=  1/(2x8192)  sec.  The  time 
histories  for  the  two  integration  time  steps  were  found  to  be  exactly  identical. 

Four  modes  are  thus  used  in  the  results  for  the  isotropic  plate  shown  in  the  following: 
The  time  histories,  PDF  and  PSD  of  maximum  deflection  and  strain  at  SPL  =90  and  120 
dB  and  AT=0.0  are  shown  in  Figures  5  and  6,  and  at  SPL=90,  100  and  120  dB  and 
ATcr/AT=2.0  are  shown  in  Figures  7  to  9,  respectively.  Table  4  gives  the  statistical 
moments  of  the  maximum  deflection  and  maximum  strain  responses.  The  skewness  and 
kurtosis  coefficients  are  defined  as 

Skewness  =  p3/a3  (43) 

Kurtosis  =  (p4/a4)-3  (44) 

where  pk  and  a  are  the  kth  central  moment  and  standard  deviation,  respectively. 
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Table  4.  Moments  of  the  Wmax/h  and  ey  for  the  14x10x0.04  in.  Isotropic  Plate 


SPL 

dB 

ATcr/AT 

Rms 

Mean 

in./in. 

Variance 

in2./in2. 

Skewness 

inVin3. 

Kurtosis 

in4./in4. 

Wmax/h 

90 

0.0 

0.1608 

1.241xl0'4 

0.0258 

0.0242 

-0.4320 

120 

0.0 

1.4039 

-2.628x1 0'3 

1.9714 

0.00634 

-0.8613 

90 

2.0 

0.8239 

-0.8163 

0.0124 

2.0074 

9.4296 

100 

2.0 

0.7470 

-0.1052 

0.5470 

-0.2480 

-1.3823 

120 

2.0 

1.5770 

1.999x1 0'3 

2.4873 

0.00272 

-0.7367 

Strain 

|iin./in. 

|ilin2./in2. 

in3. /in3. 

•  4  /•  4 

in  ./in  . 

90 

0.0 

1.324xl0'5 

-0.0176 

0.000175 

-0.0432 

-0.4050 

120 

0.0 

1.167xl0'4 

22.268 

0.01312 

0.5744 

0.4317 

90 

2.0 

4.035xl0'5 

-39.094 

0.000333 

0.22477 

0.1467 

100 

2.0 

6.6 12x1 0’5 

21.692 

0.003902 

0.26784 

-1.2727 

120 

2.0 

19. 02x1 0'5 

72.950 

0.031171 

-0.5501 

-0.1875 

At  the  low  90  dB  SPL,  the  plate  behaviors  basically  small  deflection  (Wmax/h=0.1608) 
random  vibration  dominated  by  the  fundamental  (1,1)  mode  as  shown  in  the  PSD  plots  of 
Figure  5.  The  probabilities  for  deflection  and  strain  are  both  closely  to  Gaussian.  The 
time  history  at  the  high  120  dB  SPL  in  Figure  6  is  clearly  a  large  deflection  (Wmax/h>1.0) 
nonlinear  random.  This  is  demonstrated  by  the  peaks  in  PSD  plots  that  they  are 
broadening  and  shifting  to  the  higher  frequency,  and  by  the  presence  of  a  non-zero  mean 
in-plane  strain  shown  in  the  strain  plots.  The  large  deviation  from  the  Gaussian  is  shown 
by  the  strain  PDF  and  the  larger  kurtosis  value  for  deflection  and  skewness  value  for 
strain  in  Table  4. 

At  combined  acoustic  and  thermal  loads,  the  panel  responses  show  the  three  distinct 
motions  of  the  following  (i)  small  deflection  random  vibration  about  one  of  the  two 
thermally  buckled  equilibrium  positions  in  Figure  7,  (ii)  snap-through  or  oil-canning 
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phenomenon  between  the  two  thermally  buckled  positions  in  Figure  8,  and  (iii)  large 
amplitude  nonlinear  random  vibration  covering  both  thermally  buckled  positions  shown 
in  Figure  9. 

At  low  90  dB  and  AT/ATcr=2.0,  the  time  histories  in  Figure  7  clearly  show  the  linear 
random  responses  about  one  of  the  thermally  buckled  positions  of  (Wmax/h)AT=  -0.8163. 
The  deflection  PSD  plot  shows  the  domination  of  the  fundamental  mode,  and  the  strain 
PSD  plot  shows  the  equally  contribution  from  the  third  mode.  Also,  note  that  the  general 
increase  of  the  panel  vibration  frequencies,  e.g.,  from  58  Hz  (Figure  5  at  AT=0)  to  86  Hz 
(Figure  7  at  AT/AT  Cr=2.0)  for  the  fundamental  mode.  As  the  SPL  increased  to  100  dB, 
the  time  histories  in  Figure  8  show  that  snap-through  motions  and  the  deflection  PDF  is 
non-Gaussian.  This  confirms  clearly  the  drawback  in  using  EL  approach  with  the 
Gaussian  response  assumption.  At  high  SPL  of  120  dB  in  Figure  9,  the  large  deflection 
RMS  Wmax/h  is  1.5770  which  covers  both  buckled  positions  of  (Wmax/h)AT=±0.8163. 
Nonlinearities  are  further  observed  by  the  broadening  and  shifting  of  the  peaks  in  the 
PSD  plots. 

3.4  Clamped  Composite  Plate 

Nonlinear  response  of  a  composite  plate  subjected  to  combined  acoustic  and  thermal 
loads  can  be  determined  using  the  present  modal  formulation.  Three  types  of  panel 
behavior  are  as  follows:  (i)  small  deflection  random  vibration  about  one  of  the  two 
thermally  buckled  equilibrium  positions,  (ii)  snap-through  or  oil  canning  phenomenon 
between  the  two  buckled  positions,  and  (iii)  large  nonlinear  random  vibration  covering 
both  thermally  buckled  positions,  can  be  predicted.  As  shown  in  the  second  example,  a 
clamped  rectangular  Graphite-Epoxy  plate  of  eight  layers  [0/45/-45/90]s  is  investigated. 
The  plate  is  of  15  by  12  by  0.048  in.  (38.1  by  30.5  by  0.12  cm)  and  the  material 
properties  are  the  following:  Ei=22.5  Msi  (155  GPa),  E2=1.17  (8.07),  Gi2=0.66  (4.55), 
p=0. 1468x1  O’3  lb-sec2/in.4  (1550  kg/m3),  vi2=0.22,  a,=-0.04xl0'6/°F  (-0.07x1 0'6/°C),  and 
<X2=T6.7xlO'6  (30.1x10'6).  A  proportional  damping  ratio  with  £,i=0.01  is  used.  The  in¬ 
plane  edges  are  immovable  and  the  plate  is  modeled  with  a  6  by  6  mesh  quarter  plate. 
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The  number  of  system  equations  in  structure  node  DOF  {Wb}  is  121.  The  system 
equations  are  reduced  to  the  modal  coordinates  using  the  lowest  four  symmetric  modes. 
The  maximum  deflection  and  the  maximum  strain  response  time  histories,  PDF  and  PSD 
at  AT/ATcr=2.0  and  SPL=90,  100,  and  120  dB,  are  shown  in  Figures  10  to  12, 
respectively.  The  time  histories,  PDF  and  PSD  plots  clearly  show  the  three  distinctive 
motions  of  linear  vibration  about  one  of  the  buckled  positions  at  low  90  dB  SPL,  oil¬ 
canning  jump  behavior  at  moderate  100  dB,  and  large  deflection  nonlinear  random 
vibration  at  high  SPL  of  120  dB.  It  is  also  interesting  to  note  that  the  maximum  thermal 
deflection  (Wmax/h)AT  of  the  plate  at  AT/ATcr=2.0  can  be  also  obtained  from  the  deflection 
time  histories  in  Figures  10  and  11.  The  maximum  thermal  deflection  from  the  time 
history  plots  is  Wmax/h=± 1.1446  at  AT/ATcr=2.0,  and  the  thermal  postbuckling  analyses 
are  ATcr  =36.52  °F  and  (Wmax/h)=  ±1.134  from  reference  [50].  The  deflection  PDF  plots 
shown  in  Figures  10  to  12  demonstrate  the  large  deviation  from  Gaussian  distribution. 
The  statistical  moments  for  the  maximum  deflection  and  maximum  strain  response  are 
presented  in  Table  5.  The  nonlinearities  of  the  response  are  indicated  also  by  the  PSD  in 
Figures  1 1  and  12  and  the  large  values  of  kurtosis. 


Table  5.  Moments  of  the  Wmax/h  and  £y  for  the 
Clamped  [0/45/-45/90]s  Plate  at  AT/ATCI=2.0 


SPL 

dB 

Rms 

mean 

in./in. 

variance 

in2./in2. 

skewness 

in3./in3. 

kurtosis 

in4./in4. 

90 

1.1487 

-1.1446 

Wmax/h 

0.0103 

2.4667 

18.387 

100 

1.0683 

0.4701 

0.9388 

-0.8812 

-0.8872 

120 

1.3312 

0.0064 

1.7777 

-0.0125 

-1.1648 

90 

6.107xl0'5 

|Liin./in. 

5.0390 

Strain 

|iin2./in2. 

0.003705 

in3  ./in3. 

-0.12307 

in4  ./in4. 

-0.28389 

100 

1.289xl0"4 

128.94 

0.148640 

-0.02496 

-0.26830 

120 

3. 699x1 O'4 

369.88 

1.048371 

-0.31806 

-0.10084 
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Figure  2.  Random  Noise  Generation 
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Figure  4.  Convergence  of  RMS  Maximum  Deflection  and 
Strain  for  a  Simply  Supported  Isotropic  Plate  at  120  dB  SPL 


Distribution  Range 


Distribution  Range 


Figure  5.  Response  of  a  Simply  Supported  Isotropic  Plate  at  90  dB  and  AT=0.0 


Figure  6.  Response  of  a  Simply  Supported  Isotropic  Plate  at  120  dB  and  AT=0.0 
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Figure  7.  Response  of  a  Simply  Supported  Isotropic  Plate  at  90  dB  and  AT  /ATcr=2.0 
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Figure  9.  Response  of  a  Simply  Supported  Isotropic  Plate  at  120  dB  and  AT  /ATcr=2.0 


4.  Conclusion 


A  finite  element  time  domain  modal  formulation  is  presented  for  the  prediction  of 
nonlinear  random  response  of  composite  panels  subjected  to  acoustic  pressure  at  elevated 
temperature.  The  modal  formulation  is  computationally  efficient  that  the  following 
occur:  (i)  the  number  of  modal  equations  is  small,  (ii)  the  nonlinear  modal  stiffness 
matrices  are  constant  matrices,  and  (iii)  the  time  step  of  integration  could  be  reasonably 
large.  It  is  demonstrated  that  the  following  three  types  of  panel  motions  can  be  predicted: 
(i)  linear  random  vibration  about  one  of  the  buckled  equilibrium  position,  (ii)  snap- 
through  motions  between  the  two  buckled  positions,  and  (iii)  nonlinear  random  response 
over  both  buckled  positions.  Results  of  deflection  PDF  at  the  high  SPL  also  show  that 
the  assumption  of  Gaussian  distribution  by  the  equivalent  linearization  technique  is 
inappropriate. 
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Appendix  A 


Element  Matrices  and  Load  Vectors 

Al.  Linear  Stiffness  Matrices 

[*»]=  M.Bt]dA  (Al) 

[k.Mhl  =  (A2) 

[*.]-  \jBj[A][B,]dA  (A3) 

{k„A=  lM[X*r]M<U  (A4> 

A2.  First-Order  and  Second-Order  Stiffness  Matrices 

[*u-]“  \  J,[SJKP»]<W  (A5) 

(A6) 

[*,„„]= \  far  i»n*P.]+ w 

[kn]^-[[Bem[Ale}[B,]dA  (AS) 

A3.  Load  Vectors 

0w}  =  IMiVvW  (A9> 

kl=  [[Bj{NaT}dA  (A10) 

M)}=  [{H.}p(t)dA  (All) 

A4.  Mass  Matrices 

K  ]  =  ?h\A  {Hw  lHw  J  dA  (A  12) 

kJ=  [ph{{Hu\Hv\+{HlHv\)dA  (A13) 
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Appendix  B 


Fortran  Code  for  Gaussian-Stationary  Random  Load  Generation. 

0* ********************************************************************* 


C 

C 


SIMLOAD 


£********************************************************************** 


C 

C 

C 

C 

C 

C 

C 


N 


NPT 


-NO.  OF  INTERVALS  IN  THE  SPECTRUM 
N  SHOULD  BE  AN  INTEGER  POWER  OF  TWO 
—NO. OF  POINTS  FOR  THE  TIME  SERIES 
NPT  SHOULD  BE  INTEGER  POWER  OF  TWO.  NPT>N 
I  SEED  -RANDOM  NUMBER  SEED 

TTOTAL  -  N/FMAX  TTOTAL  IS  THE  TOTAL  INTEGRATION  TIME 
DT  -  N/(NPT*FMAX)  DT  IS  THE  INTEGRATION  TIME  STEP  SIZE 


Q* ********************************************************* ************ 


C  INSTRUCTIONS  FOR  SETTING  THE  DATA 

C - 

C  1-  TAKE  HIGHEST  FREQUENCY,  FMAX 
C  2-  MINIMUM  TIME  STEP  IS  STEP_MIN=l/(2.5xFMAX) 

C  3-  N=FMAX  x  2 

C  4-  PICK  UP  TOTAL  RUNNING  TIME  (1  SEC,  2  SEC  ...)  T_total=N/FMAX 
C  5-  SELECT  NPT  TO  SATISFY  2- 
C  N 

C  STEP= - 

C  NPT  x  FMAX 

C 

0* ********************************************************************* 


PROGRAM  SIMLOAD 
IMPLICIT  REAL *8  (A-H,0-Z) 

C  COMMON  /XFER/ISTEP,DSTEP,DT,Y(16384) 

C  COMMON /XFER/DT,Y(  163  84) 

C  REAL* 8  DT,Y(2) 

DIMENSION  X(1 63  84),Y(1 6384),SP(2048),W(2048),RAND(  1 63  84) 
COMPLEX  X,ZIMAG 

OPEN  (1 , file- d:\research\load_st\pressure.dat’) 

OPEN  (2,file- d:\research\load_st\npt.dat') 

OPEN  (3,file='d:\research\load_st\fmax.dat') 

OPEN  (4  ,file- d  :\research\load_st\n  .dat') 

DATA  FMAX/1024./ 

DATA  N,NPT  /2048, 16384/ 

0* ******************************************************************** * 
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C  INITIALIZE  VARIABLES 

Q***********  *******************************************  **************** 

C  SPL=120 

C  SPP  -  8.41438*10**(-18.+SPL/10.) 

PI  =  3.1415926 
PI2  =  PI*2.0 
NP1  =N+1 

ZIMAG  =  CMPLX(0.0,1.0) 

SPPW  =  SPP/PI2 
WU  =  FMAX*PI2 
DW  =  WU/FLOAT(N) 

DO  1 19 1=1, NP1 
SP(I)  =  SPPW 
W(I)  =  (I-1)*DW 
119  CONTINUE 

AREA  =  SPP*FMAX 
SQ2DW  =  DSQRT(2.0*DW) 

TTOT  AL=PI2/D  W 
DT =TT  OT  AL/FLO  AT  (NPT) 

£********************************************************************** 
C  SET  X(1)=0.  IN  ORDER  TO  OBTAIN  NEW  MEAN  ZERO  TIME  SERIES 

^*$$ic:H:!i:$:i<:|::je:jcif:$;f:;{:$$}l::jc}{::jc:ic:i£:f:;f:;fc;i:$******************************************* 

X(  1  )=CMPLX(0 .0,0.0) 

DO  50  I=N+1,NPT 

X(I)=CMPLX(0.0,0.0) 

50  CONTINUE 

^sJj^^sjsjicsfisHsIsjIcjlcJcjfcjJcsf:}}:;};^*^^^************************************************* 

C  GENERATE  RANDOM  PHASE  ANGLES  UNIFORMLY  DISTRIBUTED 

BETWEEN  ZERO  AND  2.* PI 

£<*  sf;  3}=  ^  ;|c  s|c  sj:  3}:  sf:  s(:  3je  >}:  jJ:  *  H«  ************  *  *  *  ****  *  ***:$;*********  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  ***  *  * 

ISEED=12357 
DO  51  1=1, N 

51  RAND(I)=RAN  (I  SEED) 

DO  60  I=2,N+1 

PHI=RAND(I- 1  )*PI2 
P1=SQ2DW*DSQRT(SP(I)) 

X(I)=P1  *  CDEXP(-ZIMAG*PHI) 

60  CONTINUE 

Q*  *  *  *  *  *  *  sf;  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 

C  PERFORM  FORWARD  TRANSFORM 

Q* ************************************************************** ******* 

CALL  FFT(X,NPT,1) 

Q^^^ij}:^;^*****  *********************************************************** 

C  GET  REAL  PART 

^s|«  sf;  *****  *************************************************************  * 

DO  70 1=1, NPT 
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Y(I)=REAL(X(I)) 

70  CONTINUE 

WRITE(1,FMT=100)  Y 
100  FORMAT  (f  18. 8) 

WRITE  (2,*)  NPT 
WRITE  (3,*)  FMAX 
WRITE  (4,*)  N,DT,SPP 
STOP 
END 


Q* ^^ ********************************************************* ********* * 

C 

C  FFT 

C 

0********************************************************************** 
SUBROUTINE  FFT(X,N,K) 

IMPLICIT  INTEGER(A-Z) 

REALM  GAIN,PI2,ANG,RE,IM 
COMPLEX  X(N),XTEMP,T,U(  1 6),V,W 
LOGICAL  NEW 

DATA  PI2,GAIN,NO,KO/6.283 1 85307, 1 .0,0,0/ 


1 


2 


NEW=NO.NE.N 
IF(.N  OT.NEW)GOT  O  2 
L2N=0 
NO=l 

L2N=L2N+1 

NO=NO+NO 

IF(NO.LT.N)GOTO  1 
GAIN=1 .0/N 
AN  G=PI2  *  GAIN 
RE=COS(ANG) 

IM=SIN(ANG) 

IF(.NOT.NEW.AND.K*KO.GE.l)GOTO  4 
U  ( 1  )=CMPLX(RE,-SIGN  (IM,FLO  AT  (K))) 


DO  3  1=2, L2N 
3  U(I)=U(I-1)*U(I-1) 


KO=K 

4  SBY2=N 

DO  7  STAGE=1,L2N 
V=U(STAGE) 
W=(l. 0,0.0) 
S=SBY2 
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SBY2=S/2 
DO  6  L=1,SBY2 
DO  5  1=1, N,S 
P=I+L-1 
Q=P+SBY2 
T=X(P)+X(Q) 

X(Q)=(X(P)-X(Q))*W 

5  X(P)-T 

6  W=W*V 

7  CONTINUE 

DO  9 1=1, N 
INDEX=I-1 
JNDEX=0 
DO  8  J=1,L2N 
JNDEX=JNDEX+JNDEX 
ITEMP=INDEX/2 

IF  (ITEMP+ITEMP.NE.INDEX)  JNDEX=JNDEX+ 1 
INDEX=ITEMP 

8  CONTINUE 
J=JNDEX+1 

IF(J.LT.I)GOTO  9 
XTEMP=X(J) 

X(J)=X(I) 

X(I)=XTEMP 

9  CONTINUE 

IF(K.GT.O)RETURN 

DO  10 1=1, N 

10  X(I)=X(I)  *  GAIN 

RETURN 

END 
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Appendix  C 


Linear  Random  Vibration 

From  PDE  for  an  isotropic  rectangular  plate, 


ph^  +  DV*w  =  p0(t) 
ot 


(Cl) 


For  a  simply  supported  boundary  condition,  the  plate  deflection  and  mode  shape  are 
w(x,  y,  t)  =  £  Y,  <lmn  ('M„„  (x,y) 


f infix'^  .  ( m7ty\ 

-  x  sin  — — 

a  )  b  J 


After  substitution  of  Equation  (C2)  into  Equation  (Cl)  and  applying  the  modal 
orthogonality  condition,  the  modal  equations  are 


(C2) 


a  +n2  a  - 

Hmn  '  ™  tim'd  mn 

™,nn 


m  ,n~  1,3,5. 


Adding  a  structural  damping, 

tfnm  nm^mntf mil  + nm 

mM„ 


®mn  It 


D 

f  m ' 

2+w 

2" 

V  ph 

\aj 

•  L; 

rad/sec 


m  = 

nm 


mnK1  ph 
16 


(C3) 

(C4) 

(C5) 

(C6) 


where  co„m  and  mmn  are  the  natural  frequency  and  modal  mass,  respectively. 

The  response  to  Equation  (C4)  is  given  by  Equations  (3-57)  and  (7-37)  in  reference  [51], 
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F\n2  1-  *xSo(f) 
o  2  p  3 
tmib mn® mn 

set  mn=r  and  kl=s, 
Ek,m^ki]  =  E[qrqs]  = 


,m,  [(®r2  -  m)  J  +  4corto1  fe®,  +  Xf,®,  + 1®, ) 


The  root  mean  square  of  maximum  deflection  from  Equation  (C2)  is 


RMS(wmJ  =  \E 


M#2  ]+^.^2.  +  2.E[^i^2]h — 
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